package scu.maqiang.tiangong;

import java.util.Arrays;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES3T41;
import scu.maqiang.fes.FES3T42;
import scu.maqiang.mesh.Mesh3T4;
import scu.maqiang.mesh.MeshFormat;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.NewIterSSolver;
import scu.maqiang.numeric.SRMatrix;

public class ZKEngineHeating {

	public static void main(String[] args) {
		String dir = "TianGong\\ZKEngine\\ComputeData\\";
		String meshFile = dir + "ZKEngine.mesh";
		Mesh3T4 mesh = new Mesh3T4().readMesh(meshFile, MeshFormat.mesh);
		ScalarFunc fixedRegion = (xyz, label, param) -> xyz[0] < 0.000001?1.0:-1.0;
		int fixedLabel = 10;
		mesh.setBoundaryLabel(fixedRegion, null, fixedLabel);
		
		FES3T41 fsHeat = new FES3T41(mesh);
		int dofHeat = fsHeat.getNdof();
		double kxx = 217.0;
		double rho = 2700;
		double specificHeat = 800;
		double theta = 1.0;
		double endtime = 100;
		double dt = 5.0;
		SRMatrix A = new SRMatrix(dofHeat);
		SRMatrix M = new SRMatrix(dofHeat);
	    fsHeat.assembleStiff(new double[] {kxx}, BVPType.COMMON, A);
	    fsHeat.assembleMass(new double[] {rho * specificHeat}, BVPType.CONSIST_MASS, M);
	    
	    SRMatrix LMatrix = M.add(theta * dt, A);
		SRMatrix RMatrix = M.add((theta - 1) * dt, A); 
		
		double[] heatFlux0 = MVO.readArray(dir + "ZKEngine_Heatflux0.dat");
		double[] heatFlux10 = MVO.readArray(dir + "ZKEngine_Heatflux10.dat");
		double[] heatFluxm10 = MVO.readArray(dir + "ZKEngine_Heatflux-10.dat");
		double[] heatFlux20 = MVO.readArray(dir + "ZKEngine_Heatflux20.dat");
		double[] heatFluxm20 = MVO.readArray(dir + "ZKEngine_Heatflux-20.dat");
		
		double[] pressure0 = MVO.readArray(dir + "ZKEngine_Pressure0.dat");
		double[] pressure10 = MVO.readArray(dir + "ZKEngine_Pressure10.dat");
		double[] pressurem10 = MVO.readArray(dir + "ZKEngine_Pressure-10.dat");
		double[] pressure20 = MVO.readArray(dir + "ZKEngine_Pressure20.dat");
		double[] pressurem20 = MVO.readArray(dir + "ZKEngine_Pressure-20.dat");
		
		double E = 2.0e6;
		double Nu = 0.3;
		FES3T42 fsElastic = new FES3T42(mesh);
		int dofElastic = fsElastic.GetNdof();
		SRMatrix K = new SRMatrix(dofElastic);
		fsElastic.assembleStiff(new double[] {E, Nu}, BVPType.COMMON, K);
		
		double[] x = new double[dofHeat];
		double[] FFlux0 = new double[dofHeat];
		double time = 0.0;
		String resultFile = dir + "ZKEngineHeating.dat";
		mesh.toTecplot(resultFile, time, x);
		double[][] vecCoef = new double[1][];
		vecCoef[0] = heatFlux0;
		fsHeat.assembleFlux(vecCoef, BVPType.COMMON, FFlux0, 0, 10);
		Arrays.setAll(FFlux0, i -> FFlux0[i] * dt);
		double[] Rx = new double[dofHeat];
		double[] RHS = new double[dofHeat];
		NewIterSSolver heatSolver = new NewIterSSolver(LMatrix);
		System.out.println("Start Computing!");
		while(time < endtime) {
			time += dt;
			System.out.println("time = " + time);
			RMatrix.mul(x, Rx);
			Arrays.setAll(RHS, i -> Rx[i] + dt * FFlux0[i]);
			heatSolver.PCGSSOR(RHS, x, 1.65, 1);
			mesh.toTecplot(resultFile, time, x);
		}
	}

}
